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Abstract 

We follow up an earlier work (briefly reviewed below) to investigate the temporal stability of 
an exact travelling front solution, constructed in the form of an integral expression, for a one- 
dimensional discrete Nagumo-like model without recovery. Since the model is a piecewise linear one 
with an on-site reaction function involving a Heaviside step function, a straightforward linearisation 
around the front solution presents problems, and we follow an alternative approach in estimating 
a 'stability multiplier' by looking at the variational problem as a succession of linear evolution 
of the perturbations, punctuated with 'kicks' of small but finite duration. The perturbations get 
damped during the linear evolution, while the kicks amplify only the perturbations located at 
specific sites (the 'significant perturbations', see below) with reference to the propagating front. 
Comparison is made with results of numerical integration of the reaction-diffusion system whereby 
it appears likely that the travelling front is temporally stable for all parameter values characterising 
the model for which it exists. We modify the system by introducing a slow variation of a relevant 
recovery parameter and perform a leading order singular perturbation analysis to construct a pulse 
solution in the resulting model. In addition, we obtain (in the leading order) a 1-parameter family 
of periodic pulse trains for the system, modelling re-entrant pulses in a one-dimensional ring of 
excitable cells. 
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I. INTRODUCTION 

The discrete Nagumo (FN) model jjj without recovery on a ID lattice reads 

dur, 



dt 



D{Un+l - 2Un + Un-l) + fM, (la) 



where D is a diffusion constant coupling adjacent lattice sites and / is usually tak;en as a 
cubic bi-stable reaction function. It is considered relevant in numerous situations of interest, 
including the propagation of action potential along myelinated nerve axonSjand excitations 
of cardiac cells(see P, 0], and references in ^). Following an earlier work [J] we consider a 
piecewise linear (PWL) version of ()la|l . with the reaction function given by 

f(^u) = -u-w + eiu-a), (lb) 

where stands for the Heaviside step function, the parameter a (< 1/2) is an effective 
threshold characterising each site, separating the equilibrium values and 1 of the order 
parameter and w is a. 'recovery' variable (see below). The model (fTaj) . (fTB|) is essentially 
similar to the cubic Nagumo model in that each site is bi-stable in absence of the inter-site 
coupling. In [4I we gave an explicit construction of a travelling kink (or 'front') solution to 
this PWL version of the FN equations, of the form 

Un{t)=g{C), (2a) 

where C is a propagation variable defined as C = {xt+^) ? X being the speed of the propagating 
front. The profile function g was explicitly obtained as an integral expression 

g(Q = ap+ r [b{e)e'P' + b*{e)e-'P']e-^^'^^^-PUe, (2b) 
Jo 

where p = [(], the integer part of (, and 

ap=l-w--^Y {P>0), (3a) 
1 + 7 

= -w + -^-f-P (p<0), (3b) 
1 + 7 

111 1 

~ ~2^2D +11- i^cosO 1 - e-^ee-A'(i-'^^"^«) ' 

X{9) = -{2Dcos9-{2D + l)). (5) 
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Here the parameters 7, /i, v are defined as 

M.?^. (6b) 
. . J^, (6c) 

The speed x g^ts determined through a matching condition, for which we refer to 1^. An 
important observation is that, in contrast to a continuously distributed excitable medium, a 
travelling front on a discrete lattice gets •pinned (in this context, see 1,0]) as the threshold 
approaches a certain limiting value a (see IJl) where 

~a = -w^\\ - J\ — ^ — )]<--u;. (7) 

For a ID continuously distributed react ion- diffusion system with the same reaction funcion 
as in eq. (jlbj) a travelling kink solution would exist, for any given D, for all values of the 
threshold parameter in the range < a < 4 — w. 



The above explicit solution was obtained in 



by noting that, as the front propagates. 



there occur successive time intervals (equivalently, intervals in Q during which the system 
flla|) . ()lb|) evolves linearly, and the transition from one interval to the next is marked by 
the value of m„ crossing the threshold for some appropriate site n. The full solution is then 
obtained by an appropriate matching at these transition points between superpositions of 
the eigenmodes of the linear system. 

Fig. ^ shows the pulse profile g{C^ for arbitrarily chosen values of a, D, where one observes 
that g[C^ is piecewise continuous, with a discontinuous derivative at each lattice site owing 
to the presence of the Heaviside step function in the model (we choose the origin of time 
such that uq crosses the threshold a at t = and, accordingly, m„ crosses the threshold at 

In the first part of this paper we present a stability analysis of the travelling front (eq. 
fl2a|) . ()2b|) Vsee 0, Ql for an early work demonstrating the existence of stable travelling wave 
solutions for the Nagumo model), arriving at an estimate of a certain 'stability multiplier' 
p (see below) which is to be less than unity for the front solution to be stable, and show 



3 



numerically that conclusions arrived at on the basis of this multiplier are indeed conformed 
to in the actual time evolution of the FN system under consideration. The presence of the 
Heaviside function in the model makes the linearised evolution equation singular, involving 
delta functions at the transition points referred to above, and so a straightforward calculation 
of the linear growth rates cannot be attempted in the model. 

In section 2, we circumvent this difficulty by looking more closely at the evolution of per- 
turbations imposed on (j^af) . ()2b|). and estimate the stability multiplier, thereby arriving at 
the conclusion that the travelling front solution obtained in jj] is stable. 

In section 3 we indicate that the system (fTaj) . (jlb|l involves a certain symmetry, whereby an 
'anti-kink' solution is associated with a kink solution for any given w. We then modify the 
system (fTaj) . (fTE|) by including an equation representing the slow evolution of the recovery 
variable w, and perform a leading order singular perturbation analysis whereby a kink and 
a corresponding anti-kink solution are pieced together to yield a travelling pulse solution 
of the modified system. We also indicate how a pulse train solution involving an infinite 
periodic array of uniformly propagating pulses can be obtained in the model. 

Section 4 is devoted to concluding remarks, with brief mention of a future communication 
relating to the stability of the pulse solution. 
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FIG. 1: Kink profile function g{C,) showing discontinuity of derivative at integer values of the 
propagation variable C (corresponding to integers representing lattice sites at t = 0); parameter 
values are = 0, D = 1, a = 0.1382. 
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II. STABILITY OF THE TRAVELLING FRONT: THE STABILITY MULTIPLIER 



As mentioned above, a straightforward linearisation around the travelhng front solution 
presents technical problems involving the appearance of delta functions that make the lin- 
earised equations more singular than they actually are. Instead, if we look at the time 
evolution of a perturbation over the front solution in accordance with (eq. ()lap. ()lb|) ). we 
find that this involves a succession of intervals of linear evolution, punctuated by 'kicks' at 
the transition points wherein the evolution involves a difference between two step functions 
- one with and the other without the perturbation. For small perturbations one can then 
make a convenient approximation giving the evolution of the perturbation through the kick. 

More precisely, denoting by i]n{t) the perturbation over the travelling kink soln. Un{t) (eq. 
(j2aj), (j?B|) ). the time evolution of rjn is given by 

^^^^^ = D{7]n+i{t) - 2r]n{t) + r]n-i{t)) - r]n{t) + e{un{t) + r]n{t) - a) - e{un{t) - a). (8) 

For given n, Un(t) crosses the threshold a ai t = tn = Thus, assuming, for instance, 
rin{t) to be positive, we have 

Q{un + Vn- a) - Q{un - a) = 1 (9a) 

for 

a-7]n< Unit) < a, (9b) 

and zero otherwise. 

In the following, the time interval during which this inequality holds will be termed a 'kick'. 
As will be seen below, rjn increases monotonically during the kick and hence the kick begins 
when 

Un{t) = a-ri-, (10) 
ri~ being the value of rjn just before the kick. 

Now, for t sufficiently close to but less than tn, 

Unit) =a + iinitnW - tn) + 0(|t - tn\^), (H) 



where Un(t^) is used by recalling that is actually discontinuous at t = t„. Thus, for rjn 
sufficiently small, 

Unit) = a~r]~ (12a) 

for 

^ ^ - (12b) 

where we have assumed that Unitn) is not so small as to make the term 0(|t — t„p) relevant 
in jn}. 

Thus, the time interval during which (0(m„ + ?7„ — a) — 0('U„ — a)) differs from zero, ranges 
from (t„ — . ^".-s ) to tn- This time interval we have designated above as a 'kick'. The time 
evolution of rjnit), (ra = 0, ±1, ±2, . . .) is then one involving a series of kicks, punctuated by 
intervals where the difference of G-functions in (jSJ is zero, during which we have 

^ - Diir]n+iit) - 2(r7„(t) + ^^-lit)) + r/„(t) = 0. (13) 

For sufficiently small ?7„'s the kicks are short-lived, and of an impulsive nature. Thus, during 
the nth kick we have, to a good degree of approximation, 

while the other rjm^s remain almost unaffected during this short interval. In between kicks, 
the perturbation evolves according to (jlH|l . Equation ()14j) tells us that ?7„ indeed increases 
monotonically during the kick, and the values of rin just before and just after the nth kick 
are related as 

Vn =Vn+ Tn, (15) 

where Tn is the duration of the kick. 
Finally, using ()12b|) one gets. 



= (1 + -^)^;:, (16) 
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or, using the profile function g{() and noting tliat 

Unitn)=X9'iC)\(=0-, (17) 

we obtain, 



"-n^ + iTryJ""- '''' 

We note in passing tliat a straightforward linearisation through the replacement of {Q{un + 
Vn — a) — Q{un — a)) by 6{un — a)rin would give us 

= e^^Vn. (19) 

and would make the problem more singular then it actually is. In other words, even for 
small rjn, one has to take into account the non-linearity during the interval occupied by the 
kick, while retaining only the linear term in the duration of the kick. 

In a similar manner, if 77" happens to be negative, one finds 

Vt = . ^ 1 Vn- (20) 

X9'(0+) 

The difference between (fTH|) and (PUI) is once again a reflection of the fact that the non- 
linearity is relevant during the interval of the kick. 

In the following we present an approximate stability analysis by estimating what we have 
termed the 'stability multiplier' and our conclusions are essentially independent of whether 
we use p8|) or ()20|) for the amplifying effect of the kick on the perturbation. For the sake of 
simplicity we use below eq. (fTHj) in so far as the effect of the kick is concerned. 

We now look into eq. to see what happens to the perturbation in between the kicks. 
This is a system of linear equations and can be easily integrated for a time interval t = ^ 
from, say, the end of one kick and the begining of the next: 

rjr.it + r) = e-(^^+i)^ ^ 4_^(2Dr)r7^(t), (21) 

m 

where Ii stands for the Bessel function of order / with imaginary argument: 



Ii{u) = — I e'"'°'''e'''de. (22) 
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One can also check that evolution through ()13p leads to an over-all damping of the pertur- 
bation as expected and, for instance, obtain the following bound, 

J2M + T)\'<e-'^Y.\'^nit)\'. (23) 

n n 

Indeed, all the eigenvalues of the linear problem are real and negative, giving us the result 

We now piece together the results obtained above. Starting from time, say 0"*", just after 
the kick at the site n = (recall that a kick affects the perturbation at one site only, 
leaving unchanged the other sites), the perturbation decays through linear evolution till the 
next kick arrives at t_i = -. There follows the impulsive action of the kick, amplifying 
the perturbation at site n = — 1 by the factor on the right hand side of ()18j) . leaving the 
perturbation at the other sites unchanged. The process is repeated thereafter, the site of 
action of the kick being shifted successively by one lattice site. 

Note that, with the shift of the site of action of the kick the front itself, represented by 
fl2a|l ■ (j2b|l . moves through one lattice site during the linear evolution. Since the kick at 
n = occurs at to = when the front is also located at = 0, we conclude that each 
kick amplifies the perturbation precisely at that lattice site where the front is located at that 
instant. Since the linear evolution uniformly damps out the perturbations, we see that the 
most crucial factor in respect of the time evolution of the perturbation resides in the answer 
to the question: what happens to the perturbation at the site of location of the front as the 
latter propagates along the lattice ? Since the perturbations at the other lattice sites are 
not affected by the kicks, they are damped out. 

Thus, focussing on the perturbation at the location of the front (we call it the 'significant 
perturbation') as it gets shifted from site to site, we calculate the factor through which it 
gets amplified during the interval of one kick and the subsequent linear evolution till the 
arrival of the next kick, and call it the stability multiplier. The latter is made up of two 
factors, of which one is given by eq. (fTK|l . The other factor is to be obtained from eq. 1)2111 
as the effect of the linear evolution on the significant perturbation. 

One observes, for instance, that at the end of the time interval (of length r = -) from t^ 
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(end of kick at site n = 0) up to tZi (begining of kick at site n = —1) the perturbation at 
site n = —1 (the significant perturbation at time t_i) depends hnearly on the perturbations 
at the various sites at time t,} (i.e. the time of the previous kick). Among the latter, the 
one at site n = (i.e. the significant perturbation at time to) is of largest magnitude and 
happens to have the maximum effect on the perturbation at site n = — 1 at time t_i = = ^ 
because of its proximity. 

Thus, in order to have an order of magnitude estimate of the stability multiplier, we make 
the simplification of ignoring all perturbations on the right hand side of eq. (j?H) excepting 
the one with m = —1 which gives, for instance, 

r^_i = e-(2^+^)Vi(2Dr)r7o(0+). (24) 



Noting that rjo and r]_i are the significant perturbations at t = and t = r respectively, we 
conclude that the second factor in the stability multiplier arising due to the linear evolution 
IS e X Ji(^). 

Thus, finally we arrive at the required estimate of the stability multiplier 

/ 1 \ -(20+1) 2D, 

P = 1 + ^TTpY e^— Ji — 25a 



2tt 



1 / 1 \ -(2D + 1) / -■ 2D a 

— 1 + e^ir^ / e-''"''cosede. (25b) 

27rV X9'{^-)} ' ^ ' 



Making use of the multiplier p, we arrive at the stability criterion, 

p < 1. (26) 



In reality, eq. p6|) is not an exact criterion because (i) eq. (|T8|l is not an exact amplification 
factor due to a kick since it is valid for only one class of perturbations; (ii) it is based on 
an approach that looks at the evolution of the significant perturbation alone; and, (iii) so 
far as the effect of the linear evolution on the significant perturbation is concerned, eq. ()24|) 
overlooks perturbations at all sites excepting the significant one at the previous kick. 
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Still, one can look upon eq. ()25a|) as an effective stability criterion, that can be used as 
a convenient guide in assessing the linear stability of the propagating front solution (eq. 
(j2aj), (|2b])) of the FN system (eq. (fTa|) . (fTB|) ) under consideration. As we indicate below, 
numerical computations conform quantitatively to conclusions arrived at from eq. (j25ap . 



While accepting eq. ()25a|) as a stability criterion one, however, has to make a couple of 



qualifications re 



threshold, see j4|). Indeed, making use of results in and of the asymptotic properties of 
Bessel functions, one finds, in the pinning limit x = 0? 



ation to its validity in the limits x 



(pinning limit) and x oo (zero 



' (27a) 



2D _^ 



e X /j^( — j — ^ e X, (27b) 



i.e., in the pinning limit, we have. 



p^l. (28) 



On the other hand, in the limit x ^ oo, one similarly has, 

1 X 



X X 

i.e., once again. 



(29a) 



215+1 2D, D , 

e" — Ji(— ) ^ -, (29b) 



(29c) 



Thus, on the face of it, our theory makes no definite predictions regarding stability in the 
pinning limit as also the limit of of infinitely fast kinks. 

However, in both these limits, one notes that the amplification factor due to the kick given 
by (fT8|l goes to infinity. This, in fact, is an overestimate since for (7'(0^) 0, the second 
term in becomes vanishingly small and the third term, involving g"{Q^) is to be taken 
into account, and thus the expression for the kick duration 

= (30) 
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is an over-estimation. On the other hand, the damping factor due to hnear evolution does 
go to zero in both these hmits, and hence stabihty is recovered for both x ^ind x ^ oo. 

On the other hand, for front speeds neither too small nor too large, expression (j25aj) can be 
taken as a correct order-of-magnitude estimate. Using the front solution (j2aj), pH]) . one can 
compute p for given values of parameters a and D characterising the system. 

In fig. 12 we show the variation of the stability multiplier p with a for a set of different values 




0.05 0.1 0.15 0.2 0.25 0.3 0.35 



FIG. 2: Stability multiplier p as a function of the threshold parameter a for various different values 
of D. Fall of p near the limiting value a (x — > 0) is an artefact (see text). 

of D where we find that p < 1 for the entire range of parameter values except for a — 
(x oo) when one has p ^ 1 in accordance with ()29cp (one finds in this figure that p ^ 
in the pinning limit a — > a in apparent violation of (j2Hl); however, this is an artefact because 
the program used to evaluate ()27a|) yields an underestimation, failing to reproduce the 
exponential divergence). As already mentioned, ()29c|l are overestimations in relation to 
the actual value of p in these two limits. In other words, the travelling front solution obtained 
in our model is stable for all the parameter values for which it exists. This conclusion is 
confirmed from fig. El where we show the variation of a related multiplier p (see below) with 
a, again for a set of different values of D. We obtain p as follows: we impose a perturbation 
TjQ on the travelling front located at the site n = at time t = 0, and allow the perturbed 
system to evolve for a time t = ^, following the evolution through numerical integration. At 
the end of this interval, we look at the perturbation ri_i at site n = — 1. The ratio — is then 
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defined as p, which is thus a numerical estimate for the stability multiplier p determined 
from the actual time evolution of the system under consideration. 

A comparison of figures 121121 shows that the variations of p, p are similar in nature (excepting 
for small a, see above), indicating that our theoretical estimate gives qualitatively correct 
predictions relating to the stability of the propagating front, and the front solution is indeed 
stable in its entire range of existence. 

Fig. m presents results of numerical integration of (fTal) . ()lb|) for more general perturbations. 
We impose over the kink profile a perturbation spread over a few sites around n = at t = 
(recall that this includes what we have termed the significant perturbation) and look at the 
front profile after an appropriate time interval to see what has happened to the perturbation. 
It is found that, for all values of the parameters a, D for which the integration is performed 
(only a few representative ones among these are shown in fig. H} the perturbation dies down 
with time. 

In summary, we conclude that the stability multiplier p obtained above gives a reliable 
indication of the temporal stability of the travelling front solution and that the latter is, in 
all likelihood, stable for all values of the parameters a, D for which it exists. 




FIG. 3: Stability multiplier p obtained from numerical integration of the Nagumo system (see text) 
as a function of a for different values of D; compare with fig. (21 



12 



III. THE PROPAGATING ANTI-KINK, PULSE, AND PULSE-TRAIN SOLU- 
TIONS 



A. The anti-kink 

We begin by noting that equations (fTaf) . (fTE|) possess a symmetry w„ ^ 1 — 2a — Wn, Un — >■ 
2a result of which there is associated, with a travelhng front or kink solution, a 

travelling 'anti-kink' propagating with the same speed x- The latter is obtained from the 
former, equations ()2a|l . (j2bj) . by making the above transformation. Fig. depicts schemat- 
ically the level change in m„ (for a given lattice site n) for the kink and the corresponding 
anti-kink solution, the level change for the latter being shown in such a manner that both 
the kink and the anti-kink propagate in the same direction along the lattice. 

B. Fast and slow dynamics: the travelling pulse solution 

We now modify our system ()la|). ()lb|) in such a way that there occurs a slow change in the 
recovery variable VOfi clS cl result of which, after the level change in Un in the kink solution (for 

1.2 
1 

0.8 
0.6 
J 0-4 
0.2 


-0.2 
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FIG. 4: Numerical integration of the Nagumo system for an initial kink solution together with a 
perturbation; clockwise from top left (a) L> = l,a = 0.1382, r = 20, (b) Z) = 2,a = 0.1, r = 20, 
(c) D = 3,a = 0.1445, r = 10, (d) = 4,a = 0.2272, r = 10; in each frame, the profile obtained 
through numerical integration for a time r ('translated kink') has been compared with the profile 
obtained from the theoretically obtained solution with t = r ('expected kink'). 
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any given n, with Wn = 0) from a low to a high level (with reference to the threshold a), there 
occurs a slow rise in Wn and a corrsponding slow fall in m„ till there appears a rapid level 
change in m„ in accordance with the anti-kink solution matching the earlier kink solution, 
now with Wn = ^ — 2a. Thereafter, there occurs a slow change in Wn and u„ whereby both 
return to the resting value 0. Such alteration of fast and slow dynamics is actually observed 
in excitable media of interest, and provides the basis for FitzHugh-Nagumo dynamics (see, 
e.g., ^) as a prototype model for such systems. 

In an early paper, Conley P] gave a geometric argument based on a singular perturba- 
tion approach where he established the existence of a homoclinic orbit representing a pulse 
solution (in the context of a ID continuous excitable medium) on reducing the partial dif- 
ferential equations to a system of ODE's in terms of an appropriate propagation variable. 
Such pulse solutions in spatially continuous systems have since been extensively discussed 
in the literature (see, e.g., 1^, 11, 12,^^). In particular the question of stability of these 
pulses, obtained by piecing together a leading front and a trailing rear (respectively the 
'kink' and the 'anti-kink'), has been addressed in the context of the so-called restitution 
hypothesis 3, ll^ ■ 



In the following, we present a leading order singular-perturbation construction, along similar 
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FIG. 5: The (a) kink and (b) anti-kink profiles as function of C, (corresponding to integers repre- 
senting lattice sites at t = 0) showing level changes for given D{= 1) and a(= 0.1382); introduction 
of the slow dynamics connects up the kink and the anti-kink into a travelling pulse (see text). 
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lines, of the pulse solution in a ID discrete lattice in the above-mentioned modification of the 
system (fTaj) . ()lb|) obtained by introducing the slow dynamics involving the recovery variable 
Wn- More precisely, the system we now consider is 

du 

D{Un+l - 2Un + Un-l) + f{Un), (31a) 



dt 

dw. 



dt 



eun, (31b) 



where e is a small recovery parameter setting the time scale of the slow dynamics in terms 
of which we perform below the leading order singular perturbation calculation. 

We now have two pieces of inner solution describing a rapid level change in m„ (for w„ = 
and Wn = 1 — 2a for the kink and the anti-kink respectively) and two other pieces of outer 
solution involving the slow change of m„ as also of Wn- Fig. IHl shows schematically the 
variation of uq with t in which ABC depicts the leading edge and DE the trailing edge while 
CD and EF correspond to the slow dynamics. The points C, D, and E indicate the points of 
matching between the inner and outer solutions, which we work out below. Fig. [7| depicts 
the pulse solution in the Uq-Wq (the choice = is arbitrary ) phase plane, where the 
matching points are indicated once again, together with the points B, D' where Uq crosses 
the threshold with ~ and ~ 1 — 2a respectively. Note that the pulse rises from 
Mo = at t ^ — oo and finally recovers to mq = at t ^ +oo. Thus, from here on, we 
denote by g{() the kink solution described in section 1, with the recovery parameter w set 
at (the same notation will apply for the pulse-train solution as well). 

Following we introduce the propagation variable ( = xt ~^ ^ind represent the pulse 
solution as 

Unit) = uiO, (32a) 

Wn{t)=w{(), (32b) 

where 

limt^_±oou{C) = lim^^±^w{() = 0. (32c) 
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Let the values of the propagation variable at the matching points be denoted by Cc, Cd, and 
(e respectively, which we calculate below in the leading order in e. 



Starting from the kink solution (j2aj) . ()2b|) . for e = 0, w = 0, one obtains the leading 
correction as we switch on a small e: 



HO = - / giOdC 



(33) 



Noting that g{() ~ (here and in the following we use the symbol '~' to denote a leading 
order approximation in e) for ( away from and less than 0, while g{() ~ 1 for ( greater than 
and away from 0, we get the asymptotic expression for w{() for large (: 



HO 



X 



(34) 



With this as the inner solution for w for < C < Cc, and the corresponding inner solution 
for u as 



u{0-9{0. 
the owter solution satisfies, from eq. ()31ap . 



where, from ()31b|) 



m(c)~ 1-^(0- 



w{0 ~ 1 — e ^ 



(35) 



(36a) 



(36b) 




FIG. 6: no as a function of t for a typical pulse solution showing schematically the points of 
matching between the fast and slow dynamics; the choice n = for the lattice site is arbitrary. 
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Matching the inner and outer solutions, we get the value of (c from 



9iC< 



c 



1 - 



X 



(37) 



It is now necessary to obtain the asymptotic expression for g{() for large \(\ in order to solve 
for (c from above. This can be done by noting that, in the asymptotic region, the linear 
approximation holds for the evolution of Un{t), i.e., 

dgiO 



X- 



D{giC + 1) - 2giC) + giC - I)) - 9(0 + 9o, 



where go = 1 (resp. 0) for ( +oo (resp. ( — > — oo). Then, assuming 

^(C) -^0 ~ f^'^' (say), 

we have, 

Xlna = D{a + a'^) - {2D + 1), 
from which one can determine the exponent a. As an example, we have the results 



(7 = 7, for 0, 



and 



a = e X, for x-^ oo, 



(38) 



(39) 



(40) 



(41) 



(42) 
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FIG. 7: uq-wq (the choice n = for the lattice site is arbitrary) phase diagram for the pulse solution 
showing points of matching as also the points where uq crosses the threshold a; the origin (O) is 
the asymptotic point for the orbit for t ±00. 
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i.e., 



g{() = go + aj^^\ forx^O, (43) 

and 

g{C) = 90 + Pie'^ f^, for x^oo, (44) 
where a, {3 are constants whose values need not concern us here. 

Thus, in the following, we write g{Q = gQ + Act'*'! for |^| oo for arbitrarily specified Xi D, 
where the constant A is not relevant for our analysis, and where a is to be determined from 
the transcendental equation pUjl . 

Thus, finally, the matching at the point C gives (ref. eq. ^7\i ) 

1 + Aa^^ = 1 - f^, (45) 
X 

giving, 

Next, it is easy, from (j36b|l . to calculate to leading order in e the time or (^-interval from C 
to D, the latter being the point where w{Q ~ 1 — 2a. One thereby gets 

Cd-Cc~^M^). (47) 

As expected, the interval (Cc) after which the fast dynamics is succeeded by the slow dynam- 
ics, is small compared to the interval [Qd — Cc) during which the slow dynamics operates, in 
turn giving way to the fast level change from D to E. The interval corresponding to the latter 
is obtained as a sum of two sub-intervals : (a) the interval necessary for u{Q to drop from 
2a to a (point D' in fig. [7j), and (b) the interval for u{Q to drop from a to uc^\e ~ — 1 + 2a. 
However, in order to calculate these sub-intervals to the leading order in e, one has to again 
perform a matching between the outer and the inner solutions, as was done for the point C. 
One thereby obtains 

Ce-Cd--2—. 48 
ma 
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Finally, the slow dynamics takes over from E onwards, and the pulse returns to the level 
u = 0, w = over an infinite ^-interval. 



Summarising, we present below the travelling pulse solution obtained in the leading order 
of sungular perturbation calculation: 

-oo<C<Cc;: ^(0-^7(0, (49a) 
^(C) ~ - 9{()d(, (49b) 

X J —CO 

Cc < C < Cd : u{0 ~ e^^(^-^^\ (49c) 

w(C) ~ 1 - e~x(^-^c\ (49d) 

Cd<C<Ce: u{0^2a-g{C-^^^), (49e) 

w{C) ~ 1 - 2a, (49f) 

Cij<C<oo: m(C) ~ -(l-2a)e~xK-f^), (49g) 

w{C) ~ (1 - 2a)e"x(^-'^^), (49h) 

where Cc? Cd, Ce have been obtained above (see equations (jlHI), (jUj), (jiH|) ). Note that in 
this pulse solution the matching of u{() as also of w{() at (c, Cd, and (e is accurate only 
in the leading order in e, and so are the values of (c, Cd, Ce themselves. 

Fig. IHl shows an initial [t = 0) pulse profile constructed in accordance with eq. ()49a|l - 
fl49h|) as also the profile obtained by numerical integration of (j31ap . ()31b|) with this initial 
condition for a time lapse r (see caption). At the same time, it depicts the profile obtained 
from ()49a|) - ()49h|) by putting t = r. The agreement of the latter two indicates that ()49a|) - 
()49h|) does indeed constitute a travelling pulse solution of ()31a|) . ()31b|) . 



On the other hand, when a similar exercise is performed (fig. ^ with a pulse constructed 
as in (j49ap - ()49h|) but with matchings done at ^-values deviating from those in (|iHj) 
(we use changed values C'c = 2Cc, Cd - C'c = Cd - Cc, Ce - Cd = 2(Ci? - Cd)), one observes 
that the profile at t = r differs markedly from that obtained by numerical integration of the 
initial profile for a time lapse r, thereby confirming the validity of our singular perturbation 
construction. 
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C. The periodic pulse-train 



Simiar principles can be used to construct a one-parameter family of travelling pulse-trains 
where a pulse-train involves a periodic succession of pulses travelling along the lattice with a 
constant speed. Whereas a lone pulse involves a leading edge with a level change (from u = 
to M = 1) at w = and a trailing edge with a level change (from u = 2a to u = —1 + 2a) 
aX w = 1 — 2a, the leading edge of each pulse in a pulse train involves a level change at 
some non-zero value (say, W) of the recovery variable w. where W is to lie in the range 
< W < a — a where the upper limit is obtained by noting that a kink solution with its 
level change occurring at w = a — a gets pinned in the lattice. 

Fig. ITUl depicts schematically the variation of uq with t where successive level changes occur 
with alternating fast and slow dynamics, the points of matching between the inner and outer 
solutions being shown in a manner similar to fig. 6. Note that the leading edge of each pulse, 
instead of rising from uq = 0, joins up with the trailing edge of the previous pulse in the 
pulse train and thus, for each pulse in the train there are four matching points (C, D, E, 
F in fig. fTTH) instead of three (cf. fig. ^ for a lone pulse. Fig. ^2 shows the pulse train 



pulse profile, d=1 , chi=2.12, sigma=.67. tau=9 




FIG. 8: An initial profile obtained from tlie theoretically calculated pulse solution with t = 0, 
together with the profile ('translated pulse') obtained from the initial profile on numerical integra- 
tion of the FitzHugh-Nagumo system for a time lapse r, and the theoretically calculated profile 
drawn with t = t ('expected pulse'); D = l,e = .04, a = 0.67, a = 0.0864, r = 9; the translated 
and expected pulses coincide in the scale of the figure. 
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dynamics in the uq — wq phase plane where the closed curve is traversed repeatedly, and the 
matching points together with the points of threshold crossing are indicated. 

Applying the matching principles explained for the lone pulse, we obtain the pulse train 
solution as below (as explained above, we have a 1-parameter family of such solutions char- 
acterised by W): 

- Ci < C < Ci : uiO ~ g{0 - W, (50a) 
w{C) ~ W, (50b) 
Ci < C < Ci + C2 : u{0 ~ (1 - iy)e-^(^-^^\ (50c) 

«;(C)~l-n(C), (50d) 
Ci + C2 < C < 3Ci + C2 : u{0 ~ 2a - g{C - Ci - 2C2) + W, (50e) 

w{C) ~ 1 - 2a - W, (50f) 
3Ci + C2 < C < 3Ci + C2 + C3 : uiO ~ -(1 - 2a - W)e-i^'^-"^'-^'\ (50g) 



impropGr matching, d=1 , chi=1 .08, sigma=.55, tau=1 




-1 00 -50 50 1 00 1 50 200 250 300 

lattice site 

FIG. 9: Same as fig. |H1 but with an improper matching between the fast and slow dynamics (see 
text); D = l,€ = 0.009, a = 0.5533, a = 0.0982, r = 10 . 
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w(C) m(C), 



(50h) 



where 



and the periodicity is 



u{C + Z)=u{C). (50i) 
w{C + Z)=w{C), (50j) 



Y 1 - 

e 2a + W 

Cs ~ -^'^l ), (51c) 



^ = 4Ci + C2 + Cs- (51d) 



The speed x of the pulse-train depends on the initial level W of the recovery variable and 
is given, in the leading order of e, by the equation 




100 150 200 250 300 350 400 450 500 



FIG. 10: Mo as a function of t for a typical pulse-train solution showing schematically the matching 
points (similar to fig. EJ; the trailing edge of one pulse joins up with the leading edge of the next 
pulse, and thus there are an infinite succession of matching points. 
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(cf. eq. (26) of 



giving the speed of the kink solution for W = 0). 



Fig. ^1 shows a pulse-train profile (a) at t = and (b) at t = t (see caption) obtained 
from (j5(Ja|) - ( |50jD , together with (c) the profile obtained by numerical integration of (j31aj) . 
()31b|) for a time lapse r. The agreement between (b) and (c) indicates that the pulse-train 
solution obtained above is indeed a valid one. Additionally, fig. ^1 shows the result of a 
similar exercise but now with an improper matching between the inner and outer solutions 
using (^-intervals as in fig. El One observes here that the pulse-train profile at t = r differs 
substantially from the result of numerical integration from the initial profile. 

IV. CONCLUDING REMARKS 

In summary, we have, in this paper, established the temporal stability of the travelling 
front solution Paj) - (|Hc|l of the discrete reaction-diffusion system (fTaj) . (fTb|) by way of 
estimating the stability multiplier and have constructed the travelling pulse as also a one- 
parameter family of pulse-train solutions for a system including a slow variation of the 
recovery parameter, through a leading order singular perturbation analysis. The travelling 
pulse consists of a leading edge and a traihng edge corresponding to levels w = and 
w = 1 — 2a of the recovery variable, and also intervals of slow dynamics as explained above. 
A pulse-train, on the other hand, is made up of a periodic succession of pulses with the 

0.8 
0.7 
0.6 
0.5 

t °" 

i 0.3 
g 

" 0.2 
0.1 


^0.1 
-0.2 

-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 

fast variable u 

FIG. 11: uq-wq phase diagram for a pulse-train solution (similar to fig. Ej); the closed loop is 
traversed periodically. 
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periodic pulsetrain, d=.6, chi=.34, sigma=.36, tau=30 




600 



FIG. 12: Comparison of theoretically calculated and numerically computed pulse-train solution 
(similar to comparison in fig. Elfor a lone pulse); D = 0.6, a = 0.1395, = 0.08, e = 0.004, r = 
30, cj = 0.3661; once again, the 'translated' and 'expected' pulse trains coincide. 



improper matching, d=.6, chi=.34, sigma=.37, tau=10 




100 200 300 400 500 

lattice site 



FIG. 13: Comparison of theoretically calculated and numerically computed pulse-train solution 
(similar to comparison in fig. IHlfor a lone pulse), with an improper matching between the fast and 
slow dynamics (see text); parameters same as in fig. 1121 but with r = 10. 

leading edge of one pulse joined up with the trailing edge of the previous one. Each pulse- 
train solution belonging to the family is characterised by the parameter W, the value of 
the recovery variable at which the leading edge transition takes place (correspondingly, the 
trailing edge transition takes place at the value (1 — 2a — W) of the recovery variable). 

A periodic pulse-train is equivalent to a pulse propagating along a circular lattice, and can 
thus be used to model re-entrant pulses in a ring of excitable cells. Re-entrant waves have 
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recently been the focus of intense attention in understanding the origin of arrythmias and 



fibrillations in the cardiac tissue, and are t 
the break-up of spiral and scroll waves 



lought to constitute the mechanism underlying 



xe tnougnt to c 



In this context, we shall present, in a future communication, a stability analysis for the pulse 

I I 

solution following an approach adopted in [15] . There we consider the pulse to be travelling 
on a (sufficiently large) circular lattice (a pulse on a circular lattice is equivalent to a periodic 
pulse-train on a linear lattice; recall that a pulse-train is characterised by a parameter W) 
and reduce the problem of stability to that of obtaining the spectrum of a mapping in a 
function space. In deriving the stability criterion, we use a continuum version of (Hil), (fTB|) . 
while taking the lattice discreteness into consideration through the relation defining the 
speed X of the pulse in terms of the parameter W, a special feature of this relation being 
the pinning of the pulse at W = d — a. We shall examine as to what extent the lattice 
discreteness with its attendant feature of pinning leads to a tendency of deviation from the 
so-called restitution hypothesis, the latter being a criterion of stability expressed in terms 
of the 'diastolic interval' (DI) and the 'action potential duration' (see, e.g., |l4, 3, 2Q, 21|). 
In a separate communication we shall examine certain interesting features of pulse-train 
propagation on a discrete lattice in the presence of pacemaker oscillators, wherein the role 
of the pinning transition will be found to be especially crucial. 
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